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In modern day simulations of many-body systems much of the computational complexity is shifted to the 
identification of slowly changing molecular order parameters called collective variables (CV) or reaction co¬ 
ordinates. A vast array of enhanced sampling methods are based on the identification and biasing of these 
low-dimensional order parameters, whose fluctuations are important in driving rare events of interest. Here 
describe a new algorithm for finding optimal low-dimensional collective variables for use in enhanced sampling 
biasing methods like umbrella sampling, metadynamics and related methods, when limited prior static and 
dynamic information is known about the system, and a much larger set of candidate CVs is specified. The 
algorithm involves estimating the best combination of these candidate CVs, as quantified by a maximum 
path entropy estimate of the spectral gap for dynamics viewed as a function of that CV. Through multiple 
practical examples, we show how this post-processing procedure can lead to optimization of CV and several 
orders of magnitude improvement in the convergence of the free energy calculated through metadynamics, 
essentially giving the ability to extract useful information even from unsuccessful metadynamics runs. 


I. INTRODUCTION 

With the advent of increasingly accurate force-fields 
and powerful computers. Molecular Dynamics (MD) sim¬ 
ulations have become an ubiquitous tool for studying the 
static and dynamic properties of systems across disci¬ 
plines. However, most realistic systems of interest are 
characterized by deep, multiple free energy basins sepa¬ 
rated by high barriers. The timescales associated with 
escaping such barriers can be formidably high compared 
to what is accessible with straightforward MD even with 
the most powerful computing resources. Thus in order 
to accurately characterize such landscapes with atom¬ 
istic simulations, a large number of enhanced sampling 
schemes have become popular, starting with the pioneer¬ 
ing works of Torrie, Valleau, Bennett and otherJi^tli]. 
Many of these schemes involve probing the probabil¬ 
ity distribution along selected low-dimensional collective 
variables (CVs), either through a static pre-existing bias 
or through a bias constructed on-the-fly, that enhances 
the sampling of hard to access but important regions in 
the configuration space. 

The quality, reliability, and usefulness of the sampled 
distribution is in the end deeply dependent on the qual¬ 
ity of the chosen CV. Specifically, one key assumption 
inherent in several enhanced sampling methods is that of 
time-scale separatiorP^ for efficient and accurate sam¬ 
pling, the chosen CV should encode all the relevant slow 
dynamics in the system, and any dynamics not captured 
by the CV should be relatively fast. For most practical 
applications, there are a large number of possible CVs 
that could be chosen, and it is not at all obvious how to 
construct the best low-dimensional CV or CVs for biasing 
from these various possible options. Success in enhanced 
sampling simulations has traditionally relied on an apt 
use of physical intuition to construct such low dimen¬ 
sional CVs. Identification of good low dimensional CVs 
is in fact useful not just for enhanced sampling simula¬ 


tions such as umbrella sampling and metadynamics but 
also for distributed computing techniques like Markov 
State Models (MSMjP^, allowing one to significantly im¬ 
prove the quality and reliability of the constructed kinetic 
models. Last but not the least, having an optimal low di¬ 
mensional CV can al so hel p in the building of Brownian 
dynamics type model4^^*^. Indeed, given the importance 
of this problem, there exists a range of methods that have 
been proposed to solve ifJiSHll] 

In this communication, we report a new and com¬ 
putationally efficient algorithm for designing good low¬ 
dimensional slow CVs. We suggest that the best CV is 
one with the maximum separation of times cales between 
visible slow and hidden fast processe d^^ * ^^ ( or the max¬ 
imum spectral gap. The method is named spectral gap 
optimization of order parameters (SCOOP). Note that in 
this work henceforth we refer to the best CV in the singu¬ 
lar, without loss of any generality in the treatment. The 
notion of such a timescale separation is at the core of not 
just enhanced sampling methods but also coarse-grained, 
Multiscale and projection operator methoda^H'^. 

Our algorithm involves learning the best linear or non¬ 
linear combination of given candidate CVs, as quantified 
by a maximum path entropj®^ estimate of the spectral 
gap for the dynamics of that CV. The input to the algo¬ 
rithm is any available information about the static and 
dynamic properties of the system, accumulated through 
(i) a biased simulation performed along a sub-optimal 
trial CV, possibly (but not necessarily) complemented 
by (ii) short bursts of unbiased MD runs, or (iii) by 
knowledge of experimental observables. Any type of bi¬ 
ased simulation could be used in (i), as long as it al¬ 
lows projecting the stationary probability density esti¬ 
mate on generic CVs without having to repeat the sim¬ 
ulation. Metadynamicd^ provides this functionality in 
a straightforward manner and hence it is our method of 
choice here. Given t his in formation we use the principle 
of Maximum Calibeil^SHfll to set up an unbiased master 
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equation for the dynamics of various trial CVs. Through 
a simple post-processing optimization procedure we then 
find the CV with the maximal spectral gap of the associ¬ 
ated transfer matrix. For instance, this optimization can 
be performed through a simulated annealing approach 
that maximizes the spectral gap by performing a robust 
global search in the space of trial CVs. 

Through three practical examples, we show how our 
post-processing procedure can lead to better choices of 
CVs, and to several orders of magnitude improvement 
in the convergence of the free energy calculated through 
the popular enhanced sampling technique metadynam¬ 
ics. Furthermore, the algorithm is generally applicable 
irrespective of the number of stable basins. Our algo¬ 
rithm essentially provides the much needed ability to ex¬ 
tract useful information about relevant CVs even from 
unsuccessful metadynamics runs. In addition to use in 
free energy sampling methods, the optimized CV can 
then also be used i n other methods that provide kinetic 
rate constant^^HSl yve expect this algorithm to be of 
widespread use in designing CVs for biasing during en¬ 
hanced sampling simulations, making the process signif¬ 
icantly more automatic and far less reliant on human 
intuition. 


II. THEORY 

Let us consider a molecular system with N atoms at 
temperature T. We assume there exists a large num¬ 
ber d of available order parameters with 1 d N, 
collectively referred to as {0}, such that the dynamics 
in this d—dimensional space is Markovian. These could 
be inter-molecular distanceJ^^, torsional angles, solvation 
states, nucleus size/shap^^, bond order parameterJ^ 
etc. The identification of such order parameters poses 
another complicated problem, but as ro utinely done in 
other methods aimed at optimizing CV d^^ * ^^ *^, we as¬ 
sume such order parameters are a priori known. 

There are several available biasing techniques that can 
sample the probability distribution of the space {0}, and 
even calculate the rate constants for escape from stable 
states in this spacd^. All of these techniques are feasible 
only for a very small number of CVs whose number is 
much smaller than d - typically one to three. These are 
the order parameters whose fluctuations are deemed to 
be most important for the system or process being stud¬ 
ied, and by building a fixed or time-dependent bias of 
these CVs, one should be able to determine the true un¬ 
biased probability distribution of the full space {0}. But 
how does one decide what is an optimal low-dimensional 
subset or combination of the available order parameters? 
This dimensionality reduction is of central importance to 
methods such as umbrella sampling, metadynamics and 
others, the answer to which decides the speed of conver¬ 
gence of the biased simulation, or if it it will even ever 
converge within practically useful simulation times. 

The key idea in the current work is to perform en¬ 


hanced sampling {e.g. metadynamics) with a choice of 
trial CVs, complemented by information gathered from 
short bursts of unbiased MD simulations and experimen¬ 
tal observables when available, to iterat ively impr ove the 
CVs. The maximum Caliber frameworlPSEnEZEIl^ which 
is a dynamical generalization of the hugely popular max¬ 
imum entropy frameworlP^, provides a method for ac¬ 
complishing this. 

We start by choosing a trial CV given by /{0}, where 
/ maps the space {0} onto a lower dimensional space. 
The space along this trial CV /{0} is then discretized 
in grids labeled n. This CV could be multi-dimensional, 
with n then indexing the multidimensional grids. Let 
Pn{t) denote the instantaneous probability of the sys¬ 
tem being found in grid n. For the sake of clarity, 
we assume that / is a linear combination of {0}, i.e. 
/ = Ci0i -I- C202 + ••■ + Cd0d. The treatment developed 
here applies to non-linear combinations as well which we 
show in the examples. Then, for a fixed At, we write a 
master equation: 


Apnjt) 

At 


^rn^nmPnid'') — 

( 1 ) 


where uJnm is the rate of transition from grid n to m per 
unit time. The matrix is the entirety of all these 
rates. If the dynamics of /{0} is Markovian, then the 
matrix k of transition probabilities is given for small At 
by 

k = exp(r2At) « I -1- flAt (2) 

should not depend on the value of At used in Eq. 
This provides a self-consistency check of whether or not 
the CV so generated is Markovian. In the maximum 
Caliber approach one uses all available stationary state 
and dynamical information to construct probabilities of 
micropaths. Instead of defining the entropy as a function 
of microstate probabilities as in information theory and 
statistical thermodynamic^^, one now defines an entropy 
S' as a functional of the probabilities of micropaths, which 
is ess entia lly a path integral. For the Markovian process 
of Eq|^ 


S = -T^abPakab log kab (3) 

Path ensemble averages of time-de pende nt quantities Aab 
can now be calculated as follow^^fiiao)^ where the sub¬ 
scripts a,b denote grid indices: 

(A) = Yj(if,PakabAab (4) 

The path entropy of Eq. incremented by quantities 
accounting for constraints placed by our knowledge of 
observables {(A”j,)}, and some other constra ints s uch as 
detailed balance, is collectively called Calibeil^SESI Max¬ 
imizing the Caliber is then equivalent to being least non¬ 
committal about missing dynamic and static informa¬ 
tion, with the end result being that one obtains a re¬ 
lation between the grid-to-grid rates and the stationary 
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FIG. 1: In (a), we provide the 2-d De Leon-Berne potentiaP^^ with several candidate path CVs imposed on it. 
Black circles denote the corresponding milestone^^. See SI for further details of the CVs. In (b), the corresponding 
eigenvalues Ai and A 2 (i.e. excluding the stationary eigenvalue Aq) are shown for each of these paths. As per the 
spectral gap given here by Ai — A 2 , we identify two possible good paths marked with black circles in (b) and 
correspondingly with thicker black lines in (a). Energy is in absolute units and ksT = O.I. 


probabilities as follows: 


^ab — 


At 



(5) 


Here i runs over the number of available dynamical pieces 
of information, and pi is the Lagrange multiplier for the 
associated constraint. As a special case, consider when 
the only observable at hand is the mean number of tran¬ 
sitions in observation interval At over the entire gricP^ 
along a trial CV. In this case, the above equation takes 
a particularly simple and useful form: 




( 6 ) 


Our method then involves calculating for various trial 
CVs the spectral gap of the transition probability matrix 
k, which for a ^ 6 is kab = i^Jab^t and satisfies normal¬ 
ization T,hkab = 1- Let {A} denote the set of eigenvalues 
of k, with Ao = I > Ai > A 2 .... The spectral gap is then 
defined as As — As+i, where s is the number of barriers 
apparent from the free energy estimate projected on the 
CV at hand, that are higher than a user-defined threshold 
(typically > ksT). Estimating the Lagrange multiplier 
is computationally expensive, so a first estimate for max¬ 
imizing the spectral gap is performed using Eq. where 
the Lagrange multiplier p need not be computed. Also 
note that in the limit of small At, the matrix k will be 
diagonally dominated^, and to estimate the spectral gap 
one needs only an accurate estimate of relative local free 
energies. More static or dynamical informatiorP^f^ sim¬ 
ply introduces additional Lagrange multipliers and can 
be treated through Eq. This can be done if the inten¬ 
tion is to calculate an accurate kinetic model with correct 
estimates of the dominant eigenvalues and not just the 
spectral gap. 


We are now in a position to describe the actual algo¬ 
rithm. It comprises the following two steps in a sequential 
manner, and can be improved by iterating. 

1. Perform metadynamics along a trial CV / = Ci0i-|- 

C 2 02 + + CdQd to get a crude estimate of the 

stationary density. 

2. As post-processing, perform optimization in the 
space of mixing coefficients {ci,C 2 ...C£i} to iden¬ 
tify the CV with the maximal spectral gap. The 
reweighting functionalitjP^of metadynamics allows 
projection of free energy estimates on different CVs 
with minimal computational effort, and is used to 
calculate the k matrix through Eq. We elab¬ 
orate on the optimization procedure details in the 
next section (Illustrative Examples). 

The optimization procedure gives the best CV as the 
one with highest spectral gap, given the information at 
hand. As in any maximum entropy frameworlP^, the 
better the quality of this information, the more accurate 
will be the spectral gap. But even with very poor quality 
information, as we show in the examples, the algorithm 
still leads to signihcant improvements in the CV. Fur¬ 
thermore, whether or not the CV is Markovian can also 
be checked by repeating step 2 for different time intervals 
At of observation and determining if the spectral gap is 
independent of the value of At. 


III. ILLUSTRATIVE EXAMPLES 

A. Model 2-d landscapes: The De Leon-Berne potential 

The first illustrative example for SCOOP is a model 
2-state potential introduced by De Leon and Bern^^. To 
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FIG. 2: (a) The 5-residue peptide studied in this work. The six dihedral angles are marked, (b) The output of the 
simulated annealing algorithm run separately for different 0o values (blue circles). The starting value with the trial 
choice of CV is marked with a magenta colored star, (c) The trial (magenta) and optimized (blue) mixing 
coefficients {c} for the 6 dihedrals, (d) The spectrum of eigenvalues for dynamics projected on the trial (magenta) 
and optimized (blue) CVs. A distinct improvement can be seen in the spectral gap. 


sample this landscape at temperature ksT = 0.1, we per¬ 
form metadynamics with path CVs, a class of widely used 
CVs that can capture non-local and non-linear fluctua¬ 
tions (se^^for details). Path CVs require specification of 
a series of milestones between two points in configuration 
space, where the milestones can be described in terms of 
generic order parameters. Fluctuations in the system can 
then be enhanced in the direction along and perpendic¬ 
ular to these milestones, leading to efficient exploration 
of the space. In Fig. (a) we show the 2-d potential 
along with several possible path CVs imposed on it. We 
first perform a short trial metadynamics run biasing the 
y-coordinate. By post-processing this, we generate the 
spectral gaps for various paths using Eq. §(Fig. [T](b)). 
By comparing Fig. [^a) against Fig. [^b) , it is clear how 
the path with maximum spectral gap is the minimum en¬ 
ergy pathway passing through the saddle point. In this 
case while this result could have simply been obtained 
through Nudged Elastic Band type calculation^ - the 
point is to use this example to develop intuition for the 
method. Also note that moving around the best path 


to others that are a bit distant from it, does not lead to 
much change in the spectral gap. This is consistent with 
the observation that in several enhanced sampli ng m eth- 
ods such as metadynamics or umbrella sampling^EEl^ the 
CV need not be precisely the true reacti on coo rdinate, 
as long as it has a sufficient overlap with i-p^Ml. 

In the Supplemental Information (SI), we provide a 
similar analysis on another 2-d model potential but with 
3 states. The conclusions are similar. 


B. 5-residue peptide 

Now we move to a more complex system, which has 
also been considered as a test case for new enhanced 
sampling method^^ in order to establish their useful¬ 
ness. This is the 5-residue peptide Ace — Alas — Nme 
in vacuum (see Fig. [^(a)), where there are six possibly 
relevant dihedral torsion angles. Here we ask the ques¬ 
tion: what is the best possible I-d linear combination of 
these six dihedrals that we could bias but still maximally 
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FIG. 3: (a) and (b) show trajectories obtained from metadynamics biasing the trial CV and the optimized CV 
respectively. First 20 ns of the trajectory shown in (a) was used to generate the optimized CV for (b). A very 
pronounced improvement in the enhancement of sampling can be seen with the optimized CV. 


enhance exploration of the 6 -d space comprising all the 
dihedrals? 

In this problem, for periodicity related numerical rea¬ 
sons, we bias a reference cosine defined by cos( 6 * — 0 o), 
where 0 is one of the six dihedral angles, and dg is some 
reference value whose optimal choice we do not know a 
priori. Through our algorithm we then seek to identify: 

(a) The best choice of mixing coefficients {c} to use in 
trial CV / = ci^'^ + C 2 ^'-^ + cg ^'2 + 04'!/2 -I- 05^3 -b 
C 6 'I' 3 , where we keep the euclidean norm of {c}=l, 
and for any angle ff the prime denotes the transfor¬ 
mation d !-)• 0.5 -b cos(0 — 60 ). 

(b) The best choice of 6 * 0 , kept same for all 6 dihedrals. 

We start with the trial CV where all members of {c} 
are the same subject to euclidean norm of {c}=l, and 
an arbitrary choice of 9q = 0.75 radians is taken. A 
short metadynamics run is performed biasing this trial 
CV. See supplemental information (SI) for details of the 
metadynamics and MD parameter^^, and Fig. |^(a) for 
the metadynamics trajectory used for spectral gap opti¬ 
mization. Based on the free energy estimate generated 
from this run, a simulated annealing procedure is per¬ 
formed in the space {c} for various Og values. Starting 
from the spectral gap estimated using Eq. [^for the trial 
CV, this involves executing Metropolis moves in the {c} 
space with an attempt to find the global maxima of the 
spectral gap. In Fig. (b-d) respectively, we show how 
the spectral gap is increased by the simulated anneal¬ 
ing procedure, and the corresponding best estimate of 
{c, 6 * 0 }. The algorithm suggests the minimal role of the 
angles 'I'l, 'I'a as can be seen through their relatively 
low weight^ (Fig. (c)). The spectrum of eigenvalues 
for dynamics projected on the trial (magenta) and opti¬ 
mized (blue) CVs, along with respective spectral gaps is 
provided in Fig. [^d). Fig. [^(a-b) show the metadynam¬ 


ics trajectories for the three dihedral angles $i,$ 2 ,‘h 3 
with the trial and the optimized CVs respectively. A 
very pronounced improvement in the quality of sampling 
can be seen. Fig. |^(a-c) shows the rate of convergence of 
the error of the estimated free energjB^ with respect to 
reference values from other approache^^, through meta¬ 
dynamics runs performed with each of the trial and op- 
timi zed C Vs respectively. The error metric is the same 
as iipH, and is calculated for all points within 25 kJ 
of the global minimum in the respective 1 -d free energy. 
The behavior is robust with respect to the choice of this 
threshold value. As can be seen, the optimized CV, even 
though it was obtained on the basis of a very poorly 
converged and short (20 ns) metadynamics run, leads to 
several orders of magnitudes improvement in the rate at 
which the free energies converge. Interestingly, iterating 
the algorithm with the improved I-d CV did not lead to 
much improvement in the sampling, reflecting that the 
optimized coefficients {c} are close to the best that can 
be achieved with a 1-d CV for this problem. 

IV. CONCLUSIONS 

To conclude, we have introduced a new approach 
named SCOOP (spectral gap optimization of order pa¬ 
rameters) for improving the choice of low-dimensional 
CVs for biasing in enhanced sampling in complex sys¬ 
tems. This is accomplished through the use of maximum 
Caliber based spectral gap estimates. The algorithm is 
iterative in spirit, and attempts to learn how to improve 
CVs based on available stationary and dynamic data. 
We also provide several proof-of-concept practical exam¬ 
ples to establish the potential usefulness of the method. 
For model 2-d potentials the algorithm was shown to 
yield the minimum energy pathway. For a small peptide, 
we found very significant improvement in determining 
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FIG. 4: Errors in the 1-d free energies in kJ calculated with respect to reference free energie^^^*^ using the error 
metric fronP^. Thin and thick lines denote values using the trial and optimized CVs respectively. 


the best 1-d collective variable from six possible func¬ 
tions with no ad hoc or intuition based tuning. Future 
work will use this algorithm to treat a range of prob¬ 
lems, especially involving protein-ligand unbinding. For 
instance, the displacement of water molecules and pro¬ 
tein flexibili ty are ofte n slowly varying order parameters 
in unbindin^SmUmm^ but do we really need to bias one 
or both of these for the purpose of sampling? Another 
issue to be considered in future work is can we use these 
optimized CVs to o btain reliable dynamical information 
from metadynamicd23133i^ i ncludi ng the very important 
off-rate for ligand unbindin^^^Ml 

One central limitation of this algorithm is having to 
specify possibly a large number of order parameters that 
may be important. But for many physical problems one 
does have a sense of which order parameters could be 
at work, and this is where we expect this algorithm to 
be of tremendous use. Another obvious limitation is with 
systems devoid of a time scale separatiorP^- for example, 
in glassy systems where there is an effectively continuous 
spectrum of eigenvalues with no discernible time scale 
separation. However, the dynamics of many complex and 
real-world molecular systems does thankfully show a time 
scale separation between few relevant slow modes and 
remaining fast one^^, and we expect our algorithm to be 
of help in unraveling the thermodynamics and dynamics 
in such systems. 

ACKNOWLEDGMENTS 

We would like to thank Purushottam Dixit for help¬ 
ful discussions regarding Caliber, Omar Valsson for pro¬ 
viding system set-up and reference free energies for the 
peptide, and Jed Brown for originally suggesting a spec¬ 
tral gap approach. This work was supported by grants 
from the National Institutes of Health [NIH-GM4330] 
and the Extreme Science and Engineering Discovery En¬ 
vironment (XSEDE) [TG-MCA08X002]. 

^P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Ann 
Rev Phys Chem 53, 291 (2002). 


^G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187 (1977). 
®E. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys. 
Lett. 156, 472 (1989). 

^U. H. Hansmann and Y. Okamoto, J. Comp. Chem. 14 , 1333 
(1993). 

®Voter AF (1997) Hyperdynamics: Accelerated molecular dynam¬ 
ics of infrequent events. Phys. Rev. Lett. 78:3908-3911. 

®A. Laio and M. Parrinello, Proc Natl Acad Sci 99, 12562 (2002). 
^A. Barducci, G. Bussi, and M. Parrinello, Phys Rev Lett 100 , 
020603 (2008). 

®E. Darve, D. Rodn'guez-Gomez, and A. Pohorille, J. Chem. Phys. 
128 , 144120 (2008). 

®C. F. Abrams and E. Vanden-Eijnden, Proc. Natl. Acad. Sci. 
107 , 4961 (2010). 

^®L. Zheng, M. Chen, and W. Yang, Proc. Natl. Acad. Sci. 105 , 
20227 (2008). 

'^^P. Tiwary and A. van de Walle, Phys. Rev. B 87 , 094304 (2013) 
^^A. Berezhkovskii and A. Szabo, J. Chem. Phys. 135 , 074108 
( 2011 ). 

'^^G. Perez-Hernandez, F. Paul, T. Giorgino, G. De Fabritiis, and 
F. Noe, J. Chem. Phys. 139 , 015102 (2013). 

L. Ermak and J. McCammon, J. Chem. Phys. 69 , 1352 (1978). 
A. Morrone, J. Li, and B. J. Berne, J. Phys. Chem. B 116 , 
378 (2011). 

^®R. B. Best and G. Hummer, Proc. Natl. Acad. Sci. 102, 6732 
(2005). 

'^^R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, 
F. Warner, and S. W. Zucker, Proc. Natl. Acad. Sci. 102, 7426 
(2005). 

^®B. Peters and B. L. Trout, J. Chem. Phys. 125 , 054108 (2006). 
'^^A. Ma and A. R. Dinner, The Journal of Physical Chemistry B 
109 , 6769 (2005). 

^®M. A. Rohrdanz, W. Zheng, M. Maggioni, and C. dementi, J. 
Chem. Phys. 134, 124116 (2011). 

Ceriotti, G. A. Tribello, and M. Parrinello, Proc. Natl. Acad. 
Sci. 108 , 13023 (2011). 

Chen, T.-Q. Yu, and M. E. Tuckerman, Proc. Natl. Acad. 
Sci. 112 , 3235 (2015). 

Salvalaglio, P. Tiwary, and M. Parrinello, J. Chem. Theor. 
Comp. 10, 1420 (2014) 

'^'‘R. R. Coifman, 1. G. Kevrekidis, S. Lafon, M. Maggioni, and 
B. Nadler, Mult. Mod. Sim. 7 , 842 (2008). 

^®B. J. Berne and R. Pecora, Dynamic light scattering (Courier 
Corporation, 2000). 

^®R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). 

271. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevreki- 
did, O. Runborg, C. Theodoropoulos, Comm. Math. Sci. 1, 715 
(2003). 

2®S. Presse, K. Ghosh, J. Lee, and K. A. Dill, Rev. Mod. Phys. 
85, 1115 (2013) 





















7 


^®P. Tiwary and M. Parrinello, J Phys Chem B (2015), 
10.1021/jp504920s 

^®P. D. Dixit, A. Jain, G. Stock, and K. Dill, arXiv preprint 
arXiv:1504.01277 (2015). 

De Leon and B. Berne, J. Chem. Phys. 75, 3495 (1981). 

®®D. Branduardi, F. L. Gervasio, and M. Parrinello, J Chem Phys 
126, 054103 (2007). 

®®P. Tiwary and M. Parrinello, Phys Rev Lett 111, 230602 (2013) 
®^P. Tiwary, J. Mondal, J. A. Morrone, and B. J. Berne, Proc. 

Natl. Acad. Sci. (2015), 10.1073/pnas.1516652112 
^^P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J. Chem. 
Phys. 110, 1591 (1999). 

®®P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 
28, 784 (1983). 

T. Jaynes, Ann. Rev. Phys. Chem. 31, 579 (1980). 

®®P. D. Dixit, A. Jain, G. Stock, and K. A. Dill, J. Chem. Theor. 
Comp. (2015). 

3®E. T. Jaynes, Phys. Rev. 106, 620 (1957). 

^®A. Filyukov and V. Y. Karpov, J. Engg. Phys. Thermophys. 13, 
416 (1967). 

^^E. Rosta and G. Hummer, J. Chem. Theor. Comp. 11, 276 
(2014). 

^®D. Bicout and A. Szabo, J. Chem. Phys. 109, 2325 (1998). 

^®G. Hummer, N. Jour. Phys. 7, 34 (2005). 

^^E. Marinelli, E. Pietrucci, A. Laio, and S. Piana, PLoS Comput. 


Biol 5, el000452 (2009). 

"^^B. Berne, P. Pechukas, and G. Harp, J. Chem. Phys. 49, 3125 
(1968). 

^®D. Granata, C. Camilloni, M. Vendruscolo, and A. Laio, Proc. 
Natl. Acad. Sci. 110, 6817 (2013). 

Bonomi, C. Camilloni, A. Cavalli, and M. Vendruscolo, arXiv 
preprint arXiv: 1509.05684 (2015). 

^®G. Henkelman, B. P. Uberuaga, and H. Jonsson, J. Chem. Phys. 
113, 9901 (2000). 

^®P. Tiwary, V. Limongelli, M. Salvalaglio, and M. Parrinello, 
Proc. Natl. Acad. Sci. 112, E386 (2015). 

®®0. Valsson and M. Parrinello, Phys Rev Lett 113, 090601 (2014) 
®^0. Valsson and M. Parrinello, J. Chem. Theor. Comp. (2015). 
®®D. Branduardi, G. Bussi, and M. Parrinello, J. Chem. Theor. 
Comp. 8, 2247 (2012). 

®®G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, and 
G. Bussi, Comp. Phys. Comm. 185, 604 (2014). 

®^J. E. Ladbury, Chem. Bio. 3, 973 (1996). 

®®B. J. Berne, J. D. Weeks, and R. Zhou, Ann. Rev. Phys. Chem. 
60, 85 (2009). 

®®R. A. Copeland, D. L. Pompliano, and T. D. Meek, Nat Rev 
Drug Discov 5, 730 (2006). 

®^R. Zwanzig, Acc. Chem. Res. 23, 148 (1990). 

®®B. B. Machta, R. Chachra, M. K. Transtrum, and J. P. Sethna, 
Science 342, 604 (2013). 


